앞서 작업한 데이터와 선형회귀모형과 ARIMA 모형을 저장한 data/corona_fcst_list.rds 리스트 객체를 단변량 시계열 분석 작업을 위해서 풀어둔다.
library(tidyverse)
library(tidymodels)
library(timetk)
library(modeltime)
# 데이터 + 모형 ----
corona_fcst_list <- read_rds("data/corona_fcst_list.rds")
# 데이터 ----
full_tbl <- corona_fcst_list$data
# 모형 ----
wkfl_fit_lm <- corona_fcst_list$model$wkfl_fit_lm
wkfl_fit_arima <- corona_fcst_list$model$wkfl_fit_arima예측할 데이터와 모형개발에 활용할 데이터로 나눈 후에 모형개발에 활용할 데이터를 훈련/시험 데이터로 나눈다.
## 예측 데이터와 모형개발 데이터로 분리
forecast_tbl <- full_tbl %>%
filter(is.na(확진자))
history_tbl <- full_tbl %>%
filter(!is.na(확진자))
## 훈련/시험 데이터 분할
splits <- history_tbl %>%
time_series_split(date_var = 날짜,
assess = 30,
cumulative = TRUE)
splits<Analysis/Assess/Total>
<309/30/339>
훈련/시험 데이터로 잘 나눠졌는지 시각적으로 확인한다.
splits %>%
tk_time_series_cv_plan() %>%
plot_time_series_cv_plan(.date_var = 날짜,
.value = 확진자)tidymodels 생태계의 recipe 팩키지 recipe() 함수를 사용해서 피쳐 공학(feature engineering) 작업을 수행한다. 추후, 다양한 피처를 작업할 예정이라 대략적인 틀만 잡아둔다.
recipe_spec <- recipes::recipe(확진자 ~ 날짜, data = training(splits))
recipe_spec %>% prep() %>% juice()# A tibble: 309 x 2
날짜 확진자
<date> <int>
1 2020-01-23 0
2 2020-01-24 1
3 2020-01-25 0
4 2020-01-26 1
5 2020-01-27 1
6 2020-01-28 0
7 2020-01-29 0
8 2020-01-30 0
9 2020-01-31 7
10 2020-02-01 1
# ... with 299 more rows
단변량 예측 알고리즘으로 다음 세가지 대표적인 알고리즘을 적용시킨다.
대표적인 시계열 데이터 예측모형 ARIMA로 적합시킨다.
mdl_arima_spec <- arima_reg() %>%
set_engine("auto_arima")
wkfl_fit_arima <- workflows::workflow() %>%
add_recipe(recipe_spec) %>%
add_model(mdl_arima_spec) %>%
fit(training(splits))
wkfl_fit_arima== Workflow [trained] ==========================================================
Preprocessor: Recipe
Model: arima_reg()
-- Preprocessor ----------------------------------------------------------------
0 Recipe Steps
-- Model -----------------------------------------------------------------------
Series: outcome
ARIMA(0,1,1)
Coefficients:
ma1
-0.2198
s.e. 0.0519
sigma^2 estimated as 2729: log likelihood=-1654.94
AIC=3313.87 AICc=3313.91 BIC=3321.33
지수평활(exponential smoothing) 중 ETS(Error-Trend-Season)를 예측모형으로 적합시킨다.
mdl_ets_spec <- exp_smoothing(
error = "auto",
trend = "auto",
season = "auto",
damping = "auto"
) %>%
set_engine("ets")
wkfl_fit_ets <- workflows::workflow() %>%
add_recipe(recipe_spec) %>%
add_model(mdl_ets_spec) %>%
fit(training(splits))
wkfl_fit_ets== Workflow [trained] ==========================================================
Preprocessor: Recipe
Model: exp_smoothing()
-- Preprocessor ----------------------------------------------------------------
0 Recipe Steps
-- Model -----------------------------------------------------------------------
ETS(A,N,A)
Call:
forecast::ets(y = outcome, model = model_ets, damped = damping_ets)
Smoothing parameters:
alpha = 0.7826
gamma = 1e-04
Initial states:
l = 2.4739
s = 1.6161 1.1827 -17.5283 -11.5425 14.7329 4.8602
6.6788
sigma: 51.3469
AIC AICc BIC
4216.526 4217.265 4253.860
TBATS를 예측모형으로 적합시킨다.
mdl_tbats_spec <- seasonal_reg(
mode = "regression",
seasonal_period_1 = "auto"
) %>%
set_engine("tbats")
wkfl_fit_tbats <- workflows::workflow() %>%
add_recipe(recipe_spec) %>%
add_model(mdl_tbats_spec) %>%
fit(training(splits))
wkfl_fit_tbats== Workflow [trained] ==========================================================
Preprocessor: Recipe
Model: seasonal_reg()
-- Preprocessor ----------------------------------------------------------------
0 Recipe Steps
-- Model -----------------------------------------------------------------------
TBATS(1, {0,0}, -, {<7,2>})
Call: forecast::tbats(y = outcome, use.parallel = use.parallel)
Parameters
Alpha: 0.7812103
Gamma-1 Values: -0.002287667
Gamma-2 Values: 0.00209979
Seed States:
[,1]
[1,] 21.234291
[2,] 11.179662
[3,] -6.490442
[4,] 2.899411
[5,] -3.565564
Sigma: 50.80503
AIC: 4215.104
페이스북에서 개발한 Prophe 예측모형으로 적합시킨다.
mdl_prophet_spec <- prophet_reg(
) %>%
set_engine("prophet")
wkfl_fit_prophet <- workflows::workflow() %>%
add_recipe(recipe_spec) %>%
add_model(mdl_prophet_spec) %>%
fit(training(splits))
wkfl_fit_prophet== Workflow [trained] ==========================================================
Preprocessor: Recipe
Model: prophet_reg()
-- Preprocessor ----------------------------------------------------------------
0 Recipe Steps
-- Model -----------------------------------------------------------------------
PROPHET Model
- growth: 'linear'
- n.changepoints: 25
- changepoint.range: 0.8
- yearly.seasonality: 'auto'
- weekly.seasonality: 'auto'
- daily.seasonality: 'auto'
- seasonality.mode: 'additive'
- changepoint.prior.scale: 0.05
- seasonality.prior.scale: 10
- holidays.prior.scale: 10
- logistic_cap: NULL
- logistic_floor: NULL
- extra_regressors: 0
workflow 객체를 modeltime 팩키지 modeltime_table() 함수에 넣어 객체로 만든 후에, 추가로 지수평활법 EST 예측모형을 추가한다. modeltime_accuracy() 함수를 사용해서 모형 성능을 파악한다.
model_tbl <- modeltime_table(
wkfl_fit_arima,
wkfl_fit_ets,
wkfl_fit_tbats,
wkfl_fit_prophet
) %>%
update_model_description(.model_id = 1, "ARIMA") %>%
update_model_description(.model_id = 2, "ETS") %>%
update_model_description(.model_id = 3, "TBATS") %>%
update_model_description(.model_id = 4, "Prophet")
model_tbl %>%
modeltime::modeltime_accuracy(testing(splits))# A tibble: 4 x 9
.model_id .model_desc .type mae mape mase smape rmse rsq
<int> <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 1 ARIMA Test 275. 30.3 2.78 36.7 341. NA
2 2 ETS Test 277. 30.5 2.80 37.0 343. 0.0720
3 3 TBATS Test 277. 30.5 2.81 37.1 343. 0.0650
4 4 Prophet Test 603. 73.0 6.10 116. 647. 0.461
시험데이터를 통해 모형으로 예측한 값을 시각화한다. ETS MAE 값이 ARIMA 예측모형과 별다른 차이가 나고 있지 않지만, 신뢰구간이 생긴 것은 그나마 다행이다.
calibration_tbl <- model_tbl %>%
modeltime_calibrate(
new_data = testing(splits)
)
calibration_tbl %>%
modeltime_forecast(
new_data = testing(splits),
actual_data = history_tbl,
keep_data = TRUE,
conf_interval = 0.10
) %>%
plot_modeltime_forecast(
.legend_max_width = 60,
.legend_show = TRUE,
.conf_interval_show = TRUE,
.conf_interval_alpha = 0.20,
.conf_interval_fill = "lightblue",
.title = "코로나19 확진자 1개월 예측"
)앞선 모형은 history_tbl 을 훈련/시험 데이터에 대해 적합을 시킨 것이라 … 이제 시간을 확대하여 modeltime_refit() 함수를 사용해서 모형을 전체 데이터에 대해 다시 적합시킨다.
refit_tbl <- calibration_tbl %>%
modeltime_refit(data = history_tbl)마지막으로 앞서 구축된 모형을 바탕으로 현재까지 입수된 데이터를 바탕으로 한달 후를 예측한다.
refit_tbl %>%
modeltime_forecast(
new_data = forecast_tbl,
actual_data = history_tbl,
conf_interval = 0.3
) %>%
plot_modeltime_forecast(
.legend_max_width = 25,
.conf_interval_fill = "lightblue",
.conf_interval_alpha = 0.7,
.interactive = TRUE
)마지막 단계로 데이터와 모형을 저장하여 다음 단계로 나가기 위한 작업을 수행한다.
corona_fcst_list <- list(
# 데이터 ---
data = full_tbl,
# 예측모형 ----
model = list(
wkfl_fit_lm = wkfl_fit_lm,
wkfl_fit_arima = wkfl_fit_arima,
wkfl_fit_ets = wkfl_fit_ets,
wkfl_fit_tbats = wkfl_fit_tbats,
wkfl_fit_prophet = wkfl_fit_prophet
)
)
# 모형 + 데이터 저장 ----
corona_fcst_list %>%
write_rds("data/corona_fcst_list.rds")데이터 과학자 이광춘 저작
kwangchun.lee.7@gmail.com